Integrated Analysis of microRNA and RNA-Seq Reveals Phenolic Acid Secretion Metabolism in Continuous Cropping of Polygonatum odoratum

Polygonatum odoratum (Mill.) Druce is an essential Chinese herb, but continuous cropping (CC) often results in a serious root rot disease, reducing the yield and quality. Phenolic acids, released through plant root exudation, are typical autotoxic substances that easily cause root rot in CC. To better understand the phenolic acid biosynthesis of P. odoratum roots in response to CC, this study performed a combined microRNA (miRNA)-seq and RNA-seq analysis. The phenolic acid contents of the first cropping (FC) soil and CC soil were determined by HPLC analysis. The results showed that CC soils contained significantly higher levels of p-coumaric acid, phenylacetate, and caffeic acid than FC soil, except for cinnamic acid and sinapic acid. Transcriptome identification and miRNA sequencing revealed 15,788 differentially expressed genes (DEGs) and 142 differentially expressed miRNAs (DEMs) in roots from FC and CC plants. Among them, 28 DEGs and eight DEMs were involved in phenolic acid biosynthesis. Meanwhile, comparative transcriptome and microRNA-seq analysis demonstrated that eight miRNAs corresponding to five target DEGs related to phenolic acid synthesis were screened. Among them, ath-miR172a, ath-miR172c, novel_130, sbi-miR172f, and tcc-miR172d contributed to phenylalanine synthesis. Osa-miR528-5p and mtr-miR2673a were key miRNAs that regulate syringyl lignin biosynthesis. Nta-miR156f was closely related to the shikimate pathway. These results indicated that the key DEGs and DEMs involved in phenolic acid anabolism might play vital roles in phenolic acid secretion from roots of P. odoratum under the CC system. As a result of the study, we may have a better understanding of phenolic acid biosynthesis during CC of roots of P. odoratum.


Introduction
Polygonatum odoratum (Mill.) Druce, a traditional Chinese herb, is one of the rhizome plants in the Liliaceous family and has widespread medicinal functions. It provides multiple pharmacological effects on the blood, blood lipid, immune, endocrine, cardiovascular, and nervous systems [1,2]. In China, P. odoratum is widely distributed in 12 provinces, including the northeast, northwest, southeast, and central parts [3]. However, P. odoratum is also susceptible to replanting diseases, which result in a severe decline in the biomass and quality of underground tubers to continuous cropping (CC). Furthermore, numerous studies have reported that most medicinal plants were attacked by replanting diseases, which negatively impacts their growth and development, and they are even harvestless [4,5]. For example, ginseng yield and quality of Panax ginseng Meyer (Araliaceae) [5], Rehmannia glutinosa Libosch [6], and Radix pseudostellariae L [7] are seriously compromised under consecutive monoculture. The consecutive monoculture problem (CPM) is a major limiting Figure 1. Content of metabolites in phenolic acid biosynthesis pathway in FC and CC rhizosphere soil. ** represent significant at p ≤ 0.01 levels according to Student's t-test, respectively.

Analysis of DEGs Related to Phenolic Acid Metabolism
We performed high-throughput sequencing to determine gene expression profiles in P. odoratum roots in response to CMP. We searched for homologous sequences using BLASTX against NR, NT, KO, SwissProt, PFAM, GO, and KOG databases. The results showed that 68.337% could be assigned to a homolog in all five databases ( Figure 2A). According to the E-value distribution, we found that approximately 19.6% of the unigenes displayed very strong homology (E-value <1.0 × 10 −100 ) with available plant sequences (Figure 2B). As shown in Figure 2C, approximately 173,367 unigenes were associated with five top-hit species, including Elaeis guineensis, Phoenix dactylifera, Musa acuminate, Vitis vinifera, and Nelumbo nucifera.

Analysis of DEGs Related to Phenolic Acid Metabolism
We performed high-throughput sequencing to determine gene expression profiles in P. odoratum roots in response to CMP. We searched for homologous sequences using BLASTX against NR, NT, KO, SwissProt, PFAM, GO, and KOG databases. The results showed that 68.337% could be assigned to a homolog in all five databases (Figure 2A). According to the E-value distribution, we found that approximately 19.6% of the unigenes displayed very strong homology (E-value < 1.0 × 10 −100 ) with available plant sequences ( Figure 2B). As shown in Figure 2C, approximately 173,367 unigenes were associated with five top-hit species, including Elaeis guineensis, Phoenix dactylifera, Musa acuminate, Vitis vinifera, and Nelumbo nucifera.

Differential Expression Profiles of miRNAs and Functional Analysis during the Replanting of P. odoratum
To compare the different miRNA expression profiles in the CC and FCR of P. odoratum, differential expression analysis of the miRNAs was performed on the normalized read count for each identified miRNA. In P. odoratum roots, 253 conserved and 79 novel miRNAs were identified (Tables S2 and S3). Among them, 142 DEMs were defined by  To compare the different miRNA expression profiles in the CC and FCR of P. odoratum, differential expression analysis of the miRNAs was performed on the normalized read count for each identified miRNA. In P. odoratum roots, 253 conserved and 79 novel miRNAs were identified (Tables S2 and S3). Among them, 142 DEMs were defined by comparing the TPM expression value in CC vs FC (72 down-regulated, 70 up-regulated) P. odoratum roots (Table S4 and Figure 4A,B). Novel_6, novel_8, osa-miR167d-5p, mtr-miR319a-3p, ptc-miR396f, novel_9, ath-miR396a-5p, gma-miR396h, and zma-miR396g-3p were the most expressed DEMs in each sample ( Figure 4C).  To identify phenolic acid synthesis-associated miRNAs and mRNA interactions and understand their potential regulatory mechanisms in consecutive monoculture root tissues of P. odoratum, we performed association analyses between DEMs and target mRNAs related to phenolic acid synthesis. Among 641 target genes, five related to phenolic acid

Identification of the Key miRNA-mRNA Pairs Related to Phenolic Acid Synthesis during the Replanting of P. odoratum
To identify phenolic acid synthesis-associated miRNAs and mRNA interactions and understand their potential regulatory mechanisms in consecutive monoculture root tissues of P. odoratum, we performed association analyses between DEMs and target mRNAs related to phenolic acid synthesis. Among 641 target genes, five related to phenolic acid synthesis were screened, corresponding to eight DEMs ( Figure 5 and Tables S1 and S6). Some miRNAs were upregulated, with the expression levels of their target genes decreased in CC vs. FC root tissues of P. odoratum ( Figure 5 and Table S5). Ath-miR172a, ath-miR172c, novel_130, and tcc-miR172d which targeted ADT (Cluster-60288.94981) were up-regulated, and sbi-miR172f was down-regulated. Further, the ADT target gene had a higher expression level in CC vs FC root tissues of P. odoratum ( Figure 5), which is beneficial for phenylalanine synthesis. Moreover, osa-miR528-5p and mtr-miR2673a were down-regulated, whereas their common target gene [EC:1.11.1.7] (Cluster-60288.231604 and Cluster-60288.146690) was down-regulated, which could inhibit syringyl lignin formation [23]. Mtr-miR2673a, which targets a [EC:1.11.1.7] (Cluster-60288.201647), was down-regulated, while its other target AOC3 (Cluster-60288.117692, Cluster-60288.226755, Cluster-60288.146690) has a high expression level, coding for a critical enzyme involved in phenylacetate synthesis. Furthermore, a negative correlation was also found between nta-miR156f target to MALZ (Cluster-60288.201647), which has been related to the shikimate pathway ( Figure 5). These results indicate that these miRNAs may be involved in the phenolic acid synthesis. .146690) has a high expression level, coding for a critical enzyme involved in phenylacetate synthesis. Furthermore, a negative correlation was also found between nta-miR156f target to MALZ (Cluster-60288.201647), which has been related to the shikimate pathway ( Figure 5). These results indicate that these miRNAs may be involved in the phenolic acid synthesis.

Figure 5.
A heatmap displaying differentially expressed miRNAs (DEMs) with their target genes. FCC stands for the root samples of the first cropping, and CCR stands for root samples of consecutive cropping.

Regulatory Network of Phenolic Acid Biosynthesis
Based on RNA-seq data and the phenol propane biosynthesis method, we constructed the phenolic acid synthesis pathway of P. odoratum ( Figure 6). Combining the results of all genes and miRNAs with building a network showed a more direct relationship between miRNA-regulated gene expression, gene expression, and phenolic acid accumulation.

Regulatory Network of Phenolic Acid Biosynthesis
Based on RNA-seq data and the phenol propane biosynthesis method, we constructed the phenolic acid synthesis pathway of P. odoratum ( Figure 6). Combining the results of all genes and miRNAs with building a network showed a more direct relationship between miRNA-regulated gene expression, gene expression, and phenolic acid accumulation.
Plants 2023, 12, x FOR PEER REVIEW 8 of 17 Figure 6. Regulatory network of phenolic acid biosynthesis in the roots of P. odoratum. The red letters are up-regulation genes, the blue letters are down-regulation genes, and the green letters are miR-NAs. An arrow with a solid frame indicates only one step, while an arrow with a dotted frame arrow indicates more than one step. Black arrows indicate directions for the biosynthesis of phenolic acid.

Quantitative PCR (qPCR) Validation of the Key mRNA and miRNA Related to Phenolic Acid Synthesis
As a further validation of our RNA-seq results, we selected 19 unigenes for qPCR analysis. By utilizing shikimic acid, plants produce phenolic acids via the phenylpropanoid pathway. Hence, shikimic acid is the first step in the biosynthesis of phenolic acids in plants [24]. Among them, seven shikimic acid synthesis-related enzymes, namely  Figure 7A). The synthesis of phenolic acids in plants mainly involves the deamidation of phenylalanine. To a lesser degree, tyrosine produces cinnamic and/or p-coumaric acids. Following hydroxylation and methylation, cinnamic and p-coumaric acids form ferulic and caffeic acids [25,26]. AROK (60288.274987) and AROC (60288.232788), which encode enzymes that catalyze shikimate to phenylalanine and tyrosine, were up-regulated in CCR compared to CFR ( Figure 7A). Similarly, phenylacetate biosynthesis (AOC3 Figure 6. Regulatory network of phenolic acid biosynthesis in the roots of P. odoratum. The red letters are up-regulation genes, the blue letters are down-regulation genes, and the green letters are miRNAs. An arrow with a solid frame indicates only one step, while an arrow with a dotted frame arrow indicates more than one step. Black arrows indicate directions for the biosynthesis of phenolic acid.

Quantitative PCR (qPCR) Validation of the Key mRNA and miRNA Related to Phenolic Acid Synthesis
As a further validation of our RNA-seq results, we selected 19 unigenes for qPCR analysis. By utilizing shikimic acid, plants produce phenolic acids via the phenylpropanoid pathway. Hence, shikimic acid is the first step in the biosynthesis of phenolic acids in plants [24]. Among them, seven shikimic acid synthesis-related enzymes, namely MALZ  Figure 7A). The synthesis of phenolic acids in plants mainly involves the deamidation of phenylalanine. To a lesser degree, tyrosine produces cinnamic and/or p-coumaric acids. Following hydroxylation and methylation, cinnamic and p-coumaric acids form ferulic and caffeic acids [25,26]. AROK (60288.274987) and AROC (60288.232788), which encode enzymes that catalyze shikimate to phenylalanine and tyrosine, were up-regulated in CCR compared to CFR ( Figure 7A). Similarly, phenylacetate biosynthesis (AOC3 (60288.225005), AMIE (60288.184064)), p-coumaric acid (CYP73A, 60288.172246), caffeoyl shikimic acid  Figure 7A), whereas the lignin biosynthesis genes (E1.11.1.7, 60288.172775) decreased ( Figure 7A), which indicated that CCR was more conducive to phenolic acid synthesis, and a Pearson correlation coefficient of 0.4845 confirmed the reliability of the transcriptome sequencing results ( Figure 7C).  To validate the consistency of the miRNA and sRNA sequencing results, eight miRNAs were randomly selected for qRT-PCR analysis ( Figure 7B). The results showed that the relative expression levels of eight miRNAs correspond with those obtained from miRNA sequencing results ( Figure 7B). A significant Pearson correlation coefficient of r 2 = 0.772 indicates the reliability of by sRNA-Seq data ( Figure 7D).

Discussion
CMP, also known as replanting disease, is common to Chinese medicinal herbs in all growing regions, leading to a severe reduction in biomass and crop quality. However, the replanting disease has been attributed to many factors, including soil nutrient imbalance, physical and chemical properties, accumulation of autotoxins generated by roots, and changes in the soil microbial community structure [27,28]. Among them, the soil phenolic acids generated by roots represent an important contribution to replant disease. In addition to acidifying soil, phenolic acids induce phytotoxicity and affect soil enzyme activity, nutrient cycling, and ion uptake, thereby inhibiting the plant's growth [29,30]. Continuous cropping of strawberries, either hydroponically or with solid medium, accumulated four phenolic acids, namely p-hydroxybenzoic, ferulic, cinnamic, and p-coumaric acids [31,32]. Similarly, we previously detected five phenolic acids in the FC and CC rhizosphere soil, with significantly higher levels found in the CC soil for p-hydroxybenzoic, syringic, cumaric, and ferulic acid [13]. In the present study, we identified five other rhizosphere soil phenolic acids, three differentially accumulated in the CC and FC soil, i.e., p-coumaric acid, phenylacetate, and caffeic acid (Figure 1). These results indicated that the CC system increased phenolic acid accumulation from the rhizosphere soil of P. odoratum. Phenolic acids, as allelopathic substances, are mainly secreted by plant rhizomes but can also be decomposed by old roots and produced by soil microorganisms [18,33]. In the present study, we speculated that they might be secreted by the rhizomes of P. odoratum, thereby contributing to CMP and replant disease.
Plant phenolic acids are mainly derived via the phenylpropanoid and tyrosine pathways using shikimic acid [34]. There are three major pathways for making metabolic intermediates in plants: pentose phosphate, starch, and sucrose metabolism and pentose and glucuronate interconversions, which can provide erythrose-4-phosphate as a precursor to the synthesis of shikimic acid. In the shikimic acid biosynthesis pathway, FBP, PGD, SORD, RPE, GPI, MALZ, and GLGC are essential enzymes that produce shikimic acid [35][36][37][38]. Our previous study revealed that the pentose phosphate pathway-related genes (FBP and PGD), starch and sucrose metabolism-related gene (MALZ), and pentose and glucuronate interconversions related gene (SORD) during CC of roots of P. odoratum were significantly upregulated, which was benefic to the shikimic acid formation [13]. In this study, we found that ALDO encoded fructose-bisphosphate aldolase and was able to regulate the fructose-6-phosphate synthesis, which was used to synthesize shikimic acid ( Figure 6). Shikimic acid is the precursor of chorismate via shikimate kinase (AROK) and chorismate synthetase (AROC) in the shikimate pathway [39]. In contrast, phenolic acids are formed downstream of the chorismate, starting with the phenylpropanoid and tyrosine pathways [34]. As a precursor of phenolic acids, AROK, AROC, ADL, 4CL, CYP73A, AOC3, AMIE, and COMT are critical enzymes to produce various phenolic acids [40,41]. In this study, except for COMT, the other seven genes were markedly upregulated in the CCR ( Figure 6). As a precursor of phenolic acids, AROK encoded shikimate kinase and AROC chorismate synthase. ADT and AOC3 encoded arogenate dehydratase and primary-amine oxidase, which promote phenylacetate biosynthesis. CYP73A encoded trans-cinnamate 4-monooxygenase, which converts cinnamic acid to p-coumaric acid through hydroxylation, whereas the CoA of p-coumaric acid is esterified by 4-Couramarate CoA (4CL) [42] and then converted to caffeoyl shikimic acid or caffeic acid in the presence of shikimate O-hydroxycinnamoyl transferase (HCT) or caffeoyl shikimate esterase (CSE) (Figure 6), which is supportive of our results that p-coumaric acid, phenylacetate, and caffeic acid were significantly increased in the CC soil ( Figure 1). Conversely, the downregulation of COMT and CAD, which encoded caffeic acid 3-O-methyltransferase and cinnamyl-alcohol dehydrogenase, could inhibit lignin biosynthesis [43]. The p-coumaric acid, in turn, could be synthesized from more materials in this process. These findings further supported that the P. odoratum in the CC system might result in an accumulation of phenolic acids secreted and suggested that these genes may play essential roles in all examined phenolic acids accumulation. Another pathway for phenolic acid synthesis is mainly derived from the tyrosine pathway. Tyrosine is synthesized from chorismate by arogenate dehydrogenase (TYRAAT), which is then converted to tyramine through Tyrosine decarboxylase (TYDC) [44]. Plants can conjugate partial phenolic acids with tyramine, especially under stress [45], which activates the shikimate, phenylpropanoid, and arylmonoamine pathways [13,46], thereby accumulating more phenolic acids. Our previous study found that TYDC in CCR was significantly upregulated. In addition, we also found that TYRAAT was upregulated in this study, further suggesting that the continuous cropping of P. odoratum promoted the anabolism and accumulation of phenolic acid.
Mounting evidence suggests that miRNAs play an essential role in regulating target genes under adverse stress. We markedly identified 142 DEMs through psRobot_tar in the CC vs FC root tissues, suggesting that these miRNAs in the P. odoratum were involved in the CMP, and CMP was an extremely complex biological process. Most studies have shown that miRNAs bind to target genes' 30 untranslated regions (UTR) to suppress their expression [47]. However, in addition to sbi-miR172f overexpression down-regulating ADT, ath-miR172a, ath-miR172c, novel_130, and tcc-miR172d overexpression upregulated ADT. We inferred that after ath-miR172a, ath-miR172c, novel_130, and tcc-miR172d were overexpressed, ADT was upregulated in a compensatory manner. miRNA can act as positive regulators of target genes expression and positively regulate gene transcription by binding to targeted promoters, a process known as RNA activation [48]. This may be because a subset of th-miR172a, ath-miR172c, novel_130, and tcc-miR172d can activate ADT transcription from enhancer sites and function as an activator [49]. ADT is an arogenate dehydratase enzyme that catalyzes L-arogenate or prephenate covert to phenylalanine [50], thus increasing the phenolic acid contents of root exudates. These results implied that ath-miR172a, ath-miR172c, novel_130, tcc-miR172d, and sbi-miR172f may relate to the synthesis of phenolic acids and the increase of the CMP of the P. odoratum (Figure 6). On the contrary, osa-miR528-5p and mtr-miR2673a down-regulated their target (EC:1.11.1.7) and inhibited syringyl lignin formation [23], which could be favorable for phenolic acid accumulation. Another target of mtr-miR2673a is AOC3, which codes for a critical enzyme involved in phenylacetate synthesis, which can regulate phenolic acid synthesis in the shikimate pathway ( Figure 6). These results suggested that miRNAs may also play crucial roles in phenolic acid biosynthesis and accumulate in continuous cropping.

Plant Materials
The roots and rhizosphere soil of P. odoratum were kindly supplied by the research group of Yihong Hu [13]. The CC refers to P. odoratum grown in the land where the same plants were harvested. FC denotes P. odoratum planted the fields near CC where the cabbages had been harvested.
Following our previous method [13], rhizosphere soil was collected and preprocessed. After gentle shaking, rhizosphere soil adheres to the roots in total g. After shaking off loosely adhering soil, the rhizosphere soil was collected by brushing its circumference (0 to 5 mm), air-dried at room temperature, passed through the 2 mm mesh sieve, and stored at 4 • C until soil phenolic acid was analyzed. During the rhizome expansion stage, root samples of consecutive cropping plants (CCR) and first cropping plants (FCR) were collected randomly from the fields for microRNA and RNA-seq analysis. After collection, all samples were immediately frozen in liquid nitrogen and stored at −80 • C until further investigation.

Quantification Analysis of Phenolic Acid from Rhizosphere Soil
Rhizosphere soil phenolic acid was determined according to our previous description [13]. In brief, the rhizosphere soil was extracted using NaOH, adjusted to pH 2.5 with HCl, and extracted twice with ethyl acetate. The pooled extracts were concentrated to dryness, dissolved in methanol, and filtered through a 0.45-µm membrane. Quantification analysis of phenolic acids from rhizosphere soil was carried out using a Rigol high-performance liquid chromatography (HPLC) L300 (Rigol Technologies, Beijing, China) equipped with a Kromasil C18 column (Akzo Nobel, Amsterdam, The Netherlands). At a constant flow of 1 mL/min, 1% phosphoric acid (eluent A) and 100% methanol (eluent B) were used in a gradient. Finally, based on the retention times and adding pure standards to the sample, phenolic acids were identified and quantified at 280 nm.

RNA Isolation, Quantification, and Qualification
TRIzol reagent (Invitrogen, Carlsbad, CA, USA) was used to isolate total RNA from root tissues, which was treated with RNase-free DNase I to remove contaminants. A nanophotometer (Implen Inc., Westlake Village, CA, USA) was used to measure purity, a Qubit RNA Assay Kit (Life Technologies, Carlsbad, CA, USA) to assess concentration, and a Nano 6000 kit (Agilent Technologies, Santa Clara, CA, USA) to estimate integrity. The total RNA was immediately used for mRNA library construction, miRNA library preparation in Novogene Technology Co., Ltd. (Beijing, China), and quantitative real-time PCR (qRT-PCR) verification.

RNA-Sequencing Analysis
The FC and CC P. odoratum roots transcriptome was analyzed using two libraries (CCR and FCR) designed for RNA-Seq in our previous research [13]. Novogene Technology Co., Ltd. (Beijing, China) prepared and sequenced the library. Illumina HiSeq 2500 was used to sequence each library using paired-end protocols of 150 bp. The raw pair-end was used as clean reads.
To obtain clean reads from mRNA sequencing data, adapters, low-quality reads, and ambiguous were removed. Then, Trinity (version 2.0.6) was used for the de novo transcriptome assembly [51]. Each cluster's longest sequence region was selected as a unigene. Annotation of the assembled unigenes was performed on NR (NCBI non-redundant protein sequences database), Nt (NCBI nucleotide sequences database), Pfam (Protein family database), KOG (eukaryotic orthologous groups), COG (Clusters of Orthologous Groups), SwissProt (Swiss Institute of Bioinformatics databases), and KEGG (Kyoto Encyclopedia of Genes and Genomes databases). As a normalization method, reads per kb per million reads (RPKM) were calculated for each unigene. With the gene symbol annotation, the differential expression (n = 3) was screened using DESeq R package (Version 1.10.1) (|log2ratio| ≥ 1 and padj < 0.05) [52]. Then, the DEGs were subjected to KEGG pathway annotation using KOBAS (version 2.0). Pathways with a false discovery rate (FDR) ≤ 0.05 were significantly more enriched in DEGs [53].

Bioinformatics Analysis of Small RNA Sequences
The sequences of 18-30 nt clean reads were obtained by removing 5 adapter-contaminated reads, low-quality reads, and ambiguous reads for small RNA analysis [54]. The clean reads were then mapped against miRbase 20.0 databases (http://www.mirbase.org (10 August 2017)) [53]. The Rfam (Version 13.0) database was used to filter the rRNA, tRNA, snRNA, and snoRNA [55]. By comparing the sRNA reads with the known plant miRNAs in the miRBase 20.0, the conserved miRNAs were identified. The unannotated reads were used to predict novel miRNAs using miREvo [56].
miRNAs' expression was calculated and normalized using TPM (transcripts per kilobase million) [57]. Based on the fold difference in the expression level (|log2 fold change| ≥ 1) and the significance of the expression difference (p-value < 0.05), differentially expressed miRNAs (DEMs) in roots (FCR and CCR) were analyzed using DESeq [58].
Finally, we predicted the target genes for DEMs sequences using psRNA target [59]. Based on the corresponding genes of deferentially expressed miRNA, GO and KEGG enrichment analyses were performed with GOseq [60] and KOBAS 2.0 [57].

Quantitative Real-Time PCR Validation
To verify phenolic acid biosynthesis genes as well as miRNA, twenty genes and eight miRNAs were selected for quantitative real-time PCR (qRT-PCR) verification according to our previous methods [13]. cDNA was synthesized by a PrimerScript RT Mix (Takara, Beijing, China) using total RNA for mRNA and miRNA library construction to validate the mRNA expression level. Similarly, miRNA expression levels were validated by reversetranscribing 1 µg of the above total RNA into cDNA using miR-X miRNA First-Strand Synthesis Kit (Takara, Beijing, China). Finally, the CFX Connect Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) was used to perform the qRT-PCR analysis using the TB Green Premix Ex TaqII (Takara, Beijing, China). We use the R = 2ˆ− ∆∆Ct method to calculate each mRNA's and miRNA's relative expression levels [61]. 18S and U6 genes were used as the internal reference gene. Three technical replicates per sample were performed. Table S7 lists all primers used in this test.

Statistical Analysis
Data were expressed as the mean ± standard error of the mean (SEM). Statistical significance was assessed by the independent sample t-test. SPSS (SPSS version 21.0 for Windows; SPSS Inc., Chicago, IL, USA) software packages were used to statistical significance. There were significant differences when p < 0.05, and a tendency when 0.05 < p < 0.1. Graphs were made using GraphPad Prism version 5.01 (GraphPad Software, La Jolla, CA, USA).

Conclusions
In conclusion, in the present study, we identified five other rhizosphere soil phenolic acids, three of which were differentially accumulated in the CC and FC soil, i.e., p-coumaric acid, phenylacetate, and caffeic acid. The genomic and miRNA analysis shows for the first time the miRNA-mRNA regulatory network that may exist in phenolic acids biosynthesis. We hope that the phenolic acids biosynthesis mechanism in the continuous cropping of P. odoratum can be further improved. The results of this study will contribute to further understanding the biosynthesis of phenolic acids in the CC of roots of P. odoratum.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants12040943/s1, Figure S1: Heatmaps of phenolic acid biosynthesis pathway DEGs in the roots of P. odoratum. FCC stands for the root samples of the first cropping, and CCR stands for root samples of consecutive cropping. Table S1: Enriched DEGs relative to phenolic acid metabolism in CC vs FC root tissues of P. odoratum. Table S2: Charaterization and expression analysis of 79 novel miRNAs identified in CC vs FC root tissues of P. odoratum. Table S3: The expression of differentially expressed known miRNAs and novel miRNAs identified in CC vs FC root tissues of P. odoratum. Table S4: Target genes of known and novel miRNAs (DEMs) identified in CC vs FC root tissues of P. odoratum. Table S5: Annotation of miRNAs target genes in CC vs FC root tissues of P. odoratum. Table S6. Target genes (DEGs) of known and novel miRNAs (DEMs) involved in phenolic acid identified in CC vs FC root tissues of P. odoratum. Table S7. Primers designed for qPCR.

Conflicts of Interest:
The authors declare no conflict of interest.